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Abstract 

Two numerical procedures, one based on artificial compressibility method and the 
other pressure projection method, are outlined for obtaining time-accurate solutions of 
the incompressible Navier-Stokes equations. The performance of the two method are com- 
pared by obtaining unsteady solutions for the evolution of twin vortices behind a fiat plate. 
Calulated results are compared with experimental and other numerical results. For an un- 
steady flow which requires small physical time step, pressure projection method was found 
to be computationally efficient since it does not require any subiterations procedure. It was 
observed that the artificial compressibility method requires a fast convergence scheme at 
each physical time step in order to satisfy incompressibility condition. This was obtained 
by using a GMRES-ILU(O) solver in our computations. When a line-relaxation scheme 
was used, the time accuracy was degraded and time-accurate computations became very 
expensive. 


Introduction 

The primary objective of this research is to support the design of liquid rocket sys- 
tems for the Advanced Space Transportation System. Since the space launch systems in 
the near future are likely to rely on liquid rocket engines, increasing the efficiency and 
reliability of the engine components is an important task. One of the major problems in 
the liquid rocket engine is to understand fluid dynamics of fuel and oxidizer flows from 
the fuel tank to plume. Turbopumps in liquid rocket engines are one of the biggest source 
of vibrations. Understanding the flow through the entire turbopump geometry through 
numerical simulation will be of significant value toward design. This will help to improve 
safety of future space missions. One of the milestones of this effort is to develop, apply and 
demonstrate the capability and accuracy of 3D CFD methods as efficient design analysis 
tools on high performance computer platforms. In order to achieve fiange-to-fiange entire 
turbopump simulations, moving boundary capability and an efficient time-accurate inte- 
gration method should be build in the numerical procedure. This paper, in particular, is 
concerned with the time integration procedure of incompressible Navier-Stokes equations. 

The incompressible Navier-Stokes equations pose a special problem of satisfying the 
mass conservation equation because it is not coupled to the momentum equations. To 
satisfy incompressibility various procedures can be selected depending on the choice of for- 
mulations, variables, discretization and iterative schemes. In this paper, two formulations 
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are considered, the first one based on an artificial compressibility method and the second 
one on a pressure projection method. The artificial compressibility method 1 * takes ad- 
vantage of the advances made in conjunction with compressible flow computations. This 
approach relaxes the requirement of enforcing mass conservation equation rigorously at 
each time iteration, however, at the expense of introducing an artificial wave phenomenon. 
This approach can be viewed as a special case of a preconditioned compressible flow formu- 
lation. However, the computational efficiency is in general better than that of compressible 
flow solvers at the incompressible limit. This approach has been shown to be very robust 
in a wide range of applications 2-4 . 

The first primitive variable method for incompressible flow was developed by Harlow 
and Welch 5 using pressure projection. Numerous variants have been developed since. In 
this method, the pressure is used as a mapping parameter to satisfy the continuity equation. 
The usual computational procedure involves choosing the pressure held at the current time 
step such that continuity is satisfied at the next time step. The time step is advanced in 
multiple steps (fractional step) which is computationally convenient. However, governing 
equations are not coupled as in an artificial compressibility approach. This will affect the 
robustness and limit the maximum allowable time step size. Since this approach is time 
accurate, there are cases where the fractional step solver is computationally more efficient 
compared to the artificial compressibility method 5-9 . 

Various numerical algorithms associated with these methods have been developed 
along with accompanying flow solvers. In the present paper, it is intended to outline the 
time integration procedures of the two methods discussed above. A new time integration 
scheme is also presented for pressure projection method. Numerical results from both for- 
mulations for the development of the twin vorticies 10,11 behind the flat plate are presented 
in computed results section. 

Artificial Compressibility Formulation 

The artificial compressibility algorithm introduces a time-derivative of the pressure 
term into the continuity equation; the elliptic-parabolic type partial differential equations 
are transformed into the hyperbolic-parabolic type. The artificial compressibility method 
by Chorin (1967) can be written as 
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where t. is time, the Cartesian coordinates, u t the corresponding velocity components, 
p the pressure, ft artificial compressibility, and h, contains both convective and viscous 
terms. At steady state the pressure term in continuity equation drops out and thus in- 
compressibility is recovered. For time accurate computations, this has to be repeated at 
each time level to maintain incompressibility at each time step. 

In the present study, the time derivatives in the momentum equations are differenced 
using a second-order, three-point, backward-difference formula. 
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where u and r denote the dependent variable vector and the right hand side vector for the 
momentum equations, respectively. After the discretization in time, the pseudocompress- 
ibilty term and pseudo-time level (m) are introduced to equations. 
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Here At, A r, n, and m denote physical time step, pseudo-time step, physical time level, 
and subiteration time level, respectively. The equations are iterated to convergence in 
pseudo-time for each physical time step until the divergence of the velocity held has been 
reduced below a specified tolerance value. This typically requires 10 to 30 subiterations. 

The matrix equation is solved iteratively by using a nonfactored Gauss-Seidel type 
line-relaxation scheme, 12 which maintains stability and allows a large pseudo-time step to 
be taken. Details of the numerical method can be found in Refs. 2-3. GMRES scheme has 
also been utilized for the solution of the resulting matrix equation 13 . Computer memory 
requirement for the corresponding flow solver (INS3D-UP code) with line-relaxation is 
35 times number of grid points in words, and with GMRES-ILU(O) scheme is 220 times 
number of grid points in words. Extensive memory requirement for GMRES scheme makes 
the code unpractical for three-dimensional applications. Writing a matrix-free GMRES 
solver renmains to be one of the items for future enhancemets. 

The original version of the INS3D code 2 with pseudocompressibility approach utilized 
the Beam-Warming 14 approximate factorization algorithm and central differencing of the 
convective terms. Since the convective terms of the resulting equations are hyperbolic, 
upwind differencing can be applied to these terms. The current versions of the INS3D-UP 
code use flux-difference splitting based on the method of Roe. 15 Chakravarthy 16 outlines a 
class of high-accuracy flux-differencing schemes for the compressible flow equations. The 
third and fifth-order upwind differencing used here is an implementation of these schemes 
for the incompressible Navier-Stokes equations. The upwind differencing leads to a more 
diagonal dominant system than does central differencing and does not require a user- 
specified artificial dissipation. The viscous flux derivatives are computed by using central 
differencing. 

Time-accurate artificial compressibility formulation has been used successfully for un- 
steady calculations. The only drawback of this formulation is the computational cost due 
to subiteration procedure. 


Pressure Projection Method 

The time integration scheme is based on operator splitting, which can be accomplished 
in several ways by combining the pressure, convective, and viscous terms in the momentum 
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equations. The fractional step method is based on the decomposition of vector held into a 
divergence free component and gradient of a scalar held. The common application of this 
method is done in two steps. The first step is to solve for an auxiliary velocity held using 
the momentum equations. In the second step, the velocity held is corrected by using the 
pressure which can map the auxiliary velocity onto a divergence free velocity held. The 
momentum equations are discretized in time using a second-order implicit Runga-Kutta 
method (RK2) which can also be viewed as a predictor-corrector method. 
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where u* denotes the auxiliary velocity held. The h term in the momentum equations 
includes the convective and viscous terms. By using equation (5), equation (6) can be 
written as 
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By subtracting equation (5) from equation (7), we obtain 


^ « +1 - < ) = - Vp' + /i« +1 ) - h{u* ) (8) 

where p = p n+1 — p n . At n + 1 time level, the velocity held has to satisfy the incompress- 
ibility condition which is the continuity equation. 


V • u n+1 = 0 


( 9 ) 


This incompressiblity condition is enforced by using a Poisson equation for pressure. 
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The Poisson equation for pressure is obtained by taking the divergence of equation (8) and 
using equation (9). The only assumption is made in this procedure is that ) — h(u*) 

term in equation 8 is considered small. If the corrector step was explicit, this term would 
vanish. 

In equations (5) and (7), both convective and viscous terms are treated implicitly. 
The residual term at the (*) and (n + 1) level is linearized giving the following equations 
in delta form 
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where I is the idendity matrix. Equation (12) can also be written in more familiar form 
of fractional-step method by substituting equation (5) in equation (12). 
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Equations (11), (10), and (12) give the proposed time integration procedure of the pressure 
projection method. 

The algorithm for the pressure projection method is based on a finite-volume formu- 
lation and uses the pressure in the cell center and the mass fluxes across the faces of each 
cell as dependent variables. The discretization of the mass conservation equation in finite 
volume formulation gives 
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where S is the surface area vector. The mass conservation equation is evaluated over the 
faces of a computational Each term in equation (14) approximates the mass flux over the 
corresponding cell face. If the mass fluxes are chosen as unknowns, the continuity equation 
is satisfied automatically in generalized coordinate systems. The mass fluxes over the £, r] 
, and £ faces of the computational cell are 

U* =S s • u 

U v =S n • u (15) 

U< =S C • u 


The continuity equation with this choice of the dependent variables takes a form identical 
to the Cartesian case. Therefore, the mass fluxes are considered as the ‘natural’ depen- 
dent variables for projection methods in curvilinear coordinates. The mass conservation 
equation with new dependent variables in a generalized coordinate system becomes 
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Treating the mass fluxes as dependent variables in finite volume formulation is equivalent 
to using contravariant velocity components, scaled by the inverse of the transformation 
Jacobian, in a finite-difference formulation. This choice of mass fluxes as dependent vari- 
ables complicates the discretization of the momentum equations. In order to replace u 
by the new dependent variables U l : the corresponding area vectors are dotted with the 
momentum equations. Then the integral momentum equation is evaluated on different 
computational cells for each unknown U l . Each cell has the dimensions of x x AC, 
but the centers are located at (j + ( j,k + |,/), and (j, k, l + 4) for U^, U v , and 
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momentum equations, respectively. The staggered grid orientation eliminates pressure 
checker-board-like oscillations in pressure and provides more compact stencils. The deriva- 
tion of momentum equations and the solution procedure is outlined in reference 9. Since 
each equation is solved in a segregated fashion, memory requirements for GMRES solver 
in INS3D-FS is not as big as INS3D-UP code. Required memory for INS3D-FS is 70 times 
number of grid poins in words. 


Computed Results 

In this section, numerical results for the time evolution of twin vortices behind a flat 
plate are presented in order to verify the time-dependent features of the two algorithms. 
In order to investigate different features of the algorithms, several cases are needed to 
run with various code related parameters. To speed up this process, a two dimensional 
test case is selected here. It should be noted that associated flow solvers, INS3D-UP 
for artificial compressibility method and INS3D-FS for pressure projection method, are 
written for three-dimensional applications. With this numerical experiment, it is intended 
to give some basis for selecting a method for large three-dimensional unsteady applications 
where computing resources become a critical issue. 

Computed results from both methods are compared with the experimental data by 
Taneda and Honji 10 . The experiment has carried out in a water tank 40 cm wide. A thin 
test plate of size H = 3 cm immersed in the water was started from the rest impulsively at 
the velocity U = 0.495 cm/s. Reynolds number for this case is 126 based on U = 0.495cm/.s 
velocity and the plate height H. Computational grid with with size of 181x81x3 is presented 
in figure 1. Since INS3D-FS is written in finite volume staggered grid formulation, it 
requires one additional ghost cell in each direction. Figure 2 shows calculated velocity 
vectors obtained from INS3D-FS at various times. The flow separates the plate at each 
edge and forms a vortex pair. The twin vortices become longer in the flow direction with 
time. 

The calculated time history of the stagnation point is compared with experimental 
results and other numerical results in figure 3. Symbols represent experimental mea- 
surements, solid line and dashed line represent results from INS3D-UP and INS3D-FS, 
respectively. Dotted line show the numerical results from finite element formulations of 
Yoshida and Nomura 11 . The interval for time integration was 0.5 sec, which corresponds 
to nondimensional value of 0.0825, for all computations in figure 3. Eventhough the plate 
started impulsively in the experiment, the computations presented in figure 3 have a slow 
start procedure. Figure 4 shows prescribed velocity for an impulsive start (4a) and for 
a slow start (4b) used in INS3D-UP and INS3D-FS calculations. Reference 11 also used 
same slow start procedure in their calculations. When nondimensional time step of 0.0825 
was used with an impulsive start, large discrepancies were observed between numerical 
results and the experimental mesurements. This can be seen in figure 5a. When the time 
step is decreased, fairly good agreement was observed between numerical results and the 
measurements as seen in figure 5b. For the slow start case, the velocity profile shown in 
figure 4b is prescribed and the origin of time of calculation is appropriately shifted from 
the time of experiment. This unsteady computations with INS3D-FS (At = 0.0825) was 
completed in two hours of CPU time on single processor Cray-J90. 
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INS3D-UP computations with line-relaxation scheme is presented in figure 6. Vari- 
ous artificial compressibility parameters and number of subiterations were used. Figure 6 
shows the effects of number of subiterations and the effects of using two different artificial 
compressibility parameters (3. When the incompressibility conditions is not satisfied at 
each physical time step, numerical results can be erronous in time-accurate computations. 
With line-relaxation scheme, INS3D-UP calculations required between 4 hours of CPU 
time (10 subiterations at each physical time step) and 14 hours (40 subiterations) on a 
single processor Cray-J90 computer. Our observation for the time-accurate computations 
from this numerical example is that the artificial compressibility method requires a fast 
convergence scheme at each physical time step in order to satisfy incompressibility condi- 
tion. If this is not satisfied as seen in line-relaxation scheme, the time accuracy is degraded 
(see figure 6). In addition, artificial compressibility method with line relaxation scheme 
can be expensive for 3D time-accurate computations. In figure 7, INS3D-UP results with 
GMRES-ILU(O) solver are presented. These results were obtained less than 4 hours on a 
Cray-J90 computer. Fairly good agreement was obtained between the computed results 
and experimental data. With GMRES-ILU(O) solver, the mass flow ratio between inflow 
and exit was always satisfied. In addition, the discrepancies between numerical results 
are very small when two different values of artificial compressibility parameter were used. 
Figure 8 shows the results from artificial compressibility method with and without Poisson 
equation correction for the pressure. In artifical compressibility method, after the first 
sub-iteration, the Poisson equation is employed for the pressure correction. Chain-dashed 
line in figure 8 represents the results from this new procedure. With the Poisson equation 
correction, the line relaxation results compare well with experimental data and the GM- 
RES results with 10 subiterations. With this new procedure, both computing time and 
memory requirement are substantially reduced (at least three times). 


Concluding Remarks 

Unsteady computations were performed using two different solution algorithms, which 
are artificial compressibility method and pressure projection method. When a fast converg- 
ing scheme, such as GMRES-ILU(O) solver, was incorparated in artificial compressibility 
method, fairly good agreement was obtained between computed results and experimental 
data. Our numerical experiment showed that incompressiblity condition was satisfied in 
10 subiterations at each physical time step. Memory requirement of this scheme is the 
major drawback for three-dimensional large applications. However, memory requirement 
may not be an issue on the paralel platforms, such as SGI Origin 2000. The line-relaxation 
scheme in artificial compressibility method becomes very expensive and results in erronous 
solution for time-accurate computations. For an unsteady flow which requires small physi- 
cal time step, pressure projection method was found to be computationally efficient since it 
does not require any subiterations procedure. However, governing equations are not fully 
coupled as in the artificial compressibility approach. This may affect the robustness and 
limit the maximum allowable time step size. A new method is developed by combining 
pressure projection method with artificial compressibility method. With Poisson solver 
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correction, the number of subiteration was reduced to two iterations at each physical time 
step. 
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Figure 3 : Calculated time history 
of the stagnation point. 


Figure 1 : Computational grid for 
the flow past a 90-degree flat plate, 
(plate tickness = 0.03H) 


VELOCITY VECTORS AND MAGNITUDE CONTOURS (INS3D-FS) 


VELOCITY VECTORS AND MAGNITUDE CONTOURS (INS3D-FS) 


Figure 2a : Velocity vectors at vari- Figure 2b : Velocity vectors at vari- 
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Figure 4 : Prescribed velocity for 
an impulsive start (a) and for a slow 
start (b). 
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Figure 5a : Effects of starting pro- 
cedure. 
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Figure 5b: Effects of time-step size 
for impulsive start. 
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Figure 6 : Evaluation movement of 
stagnation point from INS3D-UP cal- 
culation with line-relaxation scheme. 
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Figure 7 : Evaluation movement of 
stagnation point from INS3D-UP cal- 
culation with GMRES-ILU(O) scheme 
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Figure 8 : Artificial compressibility 
results with and without Poisson equa- 
tion correction. 








